SHAP
install.packages(c("kernelshap", "shapviz", "ranger"))
library(kernelshap)
library(shapviz)
library(ranger)

# 2. 准备数据
X_train <- Ttrain[, !names(Ttrain) %in% "营养不良风险"]

# 3. 定义预测函数（返回概率）
p_function <- function(object, newdata) {
  predict(object, newdata = newdata, type = "prob")[, "1"]  # "1" 是阳性类别
}

# 4. 计算 SHAP 值（使用 kernelshap）
set.seed(123)
# bg_X 是背景数据集，用于计算边际均值，可随机抽取 200 行
shap_result <- permshap(
  object = rfcla,          # 您的随机森林模型
  X = X_train,             # 待解释的特征数据
  bg_X = X_train[sample(nrow(X_train), 200), ],  # 背景数据
  pred_fun = p_function    # 预测函数
)

# 5. 转换为 shapviz 对象并可视化
sv <- shapviz(shap_result)

# 5.1 特征重要性汇总图（含条形图和蜂群图）
sv_importance(sv, kind = "both", show_numbers = TRUE) +
  labs(title = "SHAP 特征重要性")
ggsave("shap_importance.png", width = 10, height = 6)

# 5.2 蜂群图
sv_importance(sv, kind = "beeswarm") +
  labs(title = "SHAP 特征贡献蜂群图")
ggsave("shap_beeswarm.png", width = 8, height = 6)

# 5.3 单个特征依赖图（如 BMI）
sv_dependence(sv, v = "BMI") +
  labs(title = "BMI 对预测的 SHAP 依赖图")
ggsave("shap_dependence_BMI.png", width = 6, height = 5)

# 5.4 单个患者解释（力图）
sv_waterfall(sv, row_id = 1) +
  labs(title = "单个患者预测的 SHAP 分解")
ggsave("shap_waterfall.png", width = 8, height = 5)

# 用您表格中的数据生成条形图的示例代码
library(ggplot2)

shap_data <- data.frame(
  特征 = c("BMI", "年龄", "三个月体重下降", "手术类型", "肠外营养", 
          "血红蛋白", "GNRI", "白蛋白", "中性粒细胞", "肠内营养"),
  SHAP值 = c(0.167, 0.137, 0.107, 0.040, 0.038, 0.020, 0.019, 0.019, 0.013, 0.005)
)

ggplot(shap_data, aes(x = reorder(特征, SHAP值), y = SHAP值, fill = SHAP值)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_gradient(low = "lightblue", high = "#D7191C") +
  labs(title = "SHAP 特征重要性",
       x = "特征", y = "平均 |SHAP 值|") +
  theme_minimal() +
  theme(legend.position = "none")
